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The synthesis of physical models for gas chemistry and turbulence from the structured 
grid codes LAURA and VULCAN into the unstructured grid code FUN3D is described. A 
directionally Symmetric, Total Variation Diminishing (STVD) algorithm and an entropy fix 
(eigenvalue limiter) keyed to local cell Reynolds number are introduced to improve solution 
quality for hypersonic aeroheating applications. A simple grid-adaptation procedure is 
incorporated within the flow solver. Simulations of flow over an ellipsoid (perfect gas, 
inviseid), Shuttle Orbiter (viscous, chemical nonequilibrium) and comparisons to the 
structured grid solvers LAURA (cylinder, Shuttle Orbiter) and VULCAN (flat plate) are 
presented to show current capabilities. The quality of heating in 3D stagnation regions is 
very sensitive to algorithm options - in general, high aspect ratio tetrahedral elements 
complicate the simulation of high Reynolds number, viscous flow as compared to locally 
structured meshes aligned with the flow. 
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speed of sound 

Jacobian of f with respect to q 
mass fraction, species s 

heat capacity per unit mass, all modes, species s 

heat capacity per unit mass, translation plus rotation modes, species s 

heat capacity per unit mass, energy mode v, species s 

laminar diffusion coefficient, species s 

turbulent diffusion coefficient, species s 

total energy per unit mass 
energy per unit mass - mode v 

average vibrational energy of preferentially dissociated species s 

inviseid flux vector 
tensile force between nodes R and L 
feature based grid adaptation function of q 
enthalpy, species s 

total enthalpy per unit mass 
thermal conductivity 
turbulent kinetic energy 
molecular weight, species s 
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degrees of freedom in translation plus rotation for species s 

unit normal and tangent vectors, respectively, to control volume face between nodes L and R 

pressure 

Prandtl number 

vector of conserved variables or convective heating 

vector of characteristic variables 

radiative energy transport term 

radius of inscribed sphere 

universal gas constant 

row and column eigenvector matrices 

distance along edge 

temperature 

time 

velocity component, j direction 
velocity components in n, £, m directions 
source term, species s 

source term, energy mode v, due to collision exchange 
Cartesian coordinate, j direction 
mole fraction, species s 

partial of pressure with respect to conserved variables pE , p r , and pe v 

Kronecker delta 

grid adaptation constant 
diagonal eigenvalue matrix 

eigenvalue of A, function of eigenvalues 

viscosity 

density 

(pcjjj) - mean stress tensor 

(px;j) - Reynolds stress tensor 

specific dissipation rate in turbulence model 


Subscripts 


i, j 

L, R 
s, r 
v, u 
t 

tr 


= x, y, or z directed momentum 
= left and right nodes 
= species 

= energy mode (vibration, electronic, ...) 
= turbulent 

= translation and rotation 
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I. Introduction 

C omputational fluid dynamics (CFD) is the numerical simulation of flowfields through the approximate solution 
of the governing partial differential equations for mass, momentum, and energy conservation coupled with the 
appropriate relations for thermodynamic and transport properties. Aerothermodynamics is the branch of fluid 
dynamics that focuses on the effects of thermodynamic and transport models on aerodynamics and heat transfer. It is 
especially focused on conditions of hypersonic velocities where the energy content and exchange between kinetic, 
internal, and chemical modes in the flow precludes the otherwise common use of calorically perfect gas 
assumptions. Computational aerothermodynamics is therefore defined in exactly the same manner as CFD, with the 
added emphasis that high temperature gas effects on pressure, skin friction, and heat transfer are included in the 
numerical simulation. The fundamental role of computational aerothermodynamics is the simulation of aerodynamic 
forces and heating for external and internal high speed flows. Reference 1 presents a review of recent applications 
for access to space and planetary missions. 

The approximate solution of the governing equations requires that the domain be subdivided into many small 
control volumes. The accuracy of the solution depends upon a variety of factors; the most critical factors are the size 
of each control volume, the orientation of its boundaries relative to a variety of flow features, and the order of 
accuracy of the discretization. At the risk of belaboring the obvious, note that hypersonic flows routinely involve 
extremes of pressure, density, and temperature separated by shocks, expansions, shear layers, and boundary layers. 
These extremes in conditions and topology of flow structure complicate computational aerothermodynamic 
simulation. The ability to orient a grid with evolving flow structures is particularly challenging. For example, 
accurate simulation of heat transfer requires adequate resolution of the boundary layer (or merged layer) and 
accurate representation of conditions at the boundary-layer edge. Conditions at the boundary-layer edge in turn, par- 
ticularly in the stagnation region of hypersonic flows, are dependent on entropy carried along streamlines from the 
shock. Any irregularities in the captured shock create associated irregularities in entropy that feed the rest of the 
domain. Further complicating the simulation is the fact that the natural stability properties of upwind schemes so 
often used in hypersonic applications are in turn most poorly realized in the broad, stagnation region of blunt bodies 
where system eigenvalues are small. 

Nevertheless, computational aerothermodynamic tools that generally do a good job in handling the difficulties of 
steady, laminar, hypersonic, blunt body and attached flow simulation have evolved over the past 15 years. The 
greatest uncertainty in such applications is the validation of physical models for a host of non-equilibrium 
processes 2 . Greater challenges are encountered in the simulation of separated flows, complex shock-shock and shock 
- boundary-layer interactions, prediction of transition, time-dependent flows, and plasma flows. The application of 
unstructured grid methods for such three-dimensional simulations (even the relatively simple blunt body problem) 
has proven to be especially challenging when heat transfer is required as will be discussed herein. Algorithms 
observed to improve simulation quality are presented, though significant problems remain. 

The synthesis of the physical models within the structured codes LAURA 3 (focus on external, hypersonic flow 
simulations) and VULCAN 4 (focus on internal, scramjet flow simulations) into the unstructured flow solver 
FUN3D 5 (perfect gas flow simulations with adjoint-based design and optimization capabilities) has proceeded 
within the FTigh Energy Flow Solver Synthesis 6 (F1EFSS) project. The ultimate goal of this project is to develop 
robust schemes with quantified uncertainties requiring minimal user intervention for aerothermodynamic 
applications. This capability exists as a generic gas model option utilizing FIEFSS modules within the code FUN3D. 


II. Grid 

The FIEFSS team is investing in unstructured grid technology, including the use of mixed elements - tetrahedron 
(tets), prisms, pyramids, and hexahedron (hexes) - as the underlying support for future hypersonic flow simulation. 
We defer any investment in patched or overset multiblock structured grid systems at this time. The rationale for this 
approach is based on our assessment 7 of current and anticipated future capabilities for grid generation and grid 
adaptation. This assessment considers: (1) efficiency in the initial generation of a “reasonable” grid over complex 
configurations; (2) ability to adapt the grid to flow structures (steady and time-dependent) and changes in 
parameterized structures on the vehicles; and (3) solution quality and efficiency on a “good” grid using best 
available algorithms. A “reasonable” grid is one that will accommodate a converged solution with at least 
approximate resolution of flow structures. It provides an initial condition for future adaptation and refinement. A 
“good" grid is one that will produce a grid converged solution such that an additional refinement will result in 
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dependent variable changes within user specified tolerances. The three criteria defined above are discussed in more 
detail in a related paper 7 - especially as regards simulations of hypersonic flows. 

Only grids composed entirely of tets are used in the present work. We forego the usual practice of resolving the 
near wall boundary layer regions with hexes or prisms. It is not clear if structured or semi-structured grids are 
required to accurately simulate convective heating and shear. (We have been unable to find published reports with a 
focus on heating or shear that have used tets all the way to the wall.) We start here with a view that a requirement 
for a semi-structured grid across boundary layers and free shear layers is a constraint that, at a minimum, 
complicates grid generation and adaptation algorithms. If we find that tets are unacceptable for boundary layer 
resolution we must invest in grid adaptation algorithms that will accommodate resolution of free shear layers (and 
possibly shock waves) with high-aspect-ratio, semi-structured cells properly aligned with the flow features. For now, 
we look to extract as much capability as possible using only tets. 

III. Physical Models 


The compressible Reynolds averaged equations for species continuity, momentum, total energy, and modal 
(vibration, electronic) energy conservation are written as: 
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is the mean stress tensor. Repeated indices indicate summation 


over direction (j or k), species (s), or energy mode (v). The Reynolds stress tensor, pz t - , is a product of the 
averaging process and must be modeled to close the system of equations. 


A. Thermodynamic Models 

Three options for the equation of state: (1) calorically perfect, (2) chemical equilibrium, and (3) chemical 
nonequilibrium mixture of thermally perfect gases are accommodated within the HEFSS thermodynamics module. 
In the case of chemical equilibrium, p(p, C) is available from curve fits for air. (Other gas mixtures with constant 
elemental mass fraction may be accommodated in like manner.) In the case of chemical nonequilibrium, a mixture 
of thermally perfect gases with p = p s RT / M s , is assumed. Curve fits for heat capacity of species s, C p s (T) , 
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and related thermodynamic properties are available in the literature. 8 ' 9 In the event that a temperature falls outside 
the curve fit range it is assumed that C p s (T) = C p s (T min ) if T < T min and C ps (T) = C p s (T max ) if 

T > T max . 

Thermal nonequilibrium is simulated with a two-temperature model 10 in which heavy particle translation and 
rotation energy modes are assumed in equilibrium at temperature T and vibration, electronic, and free electron 
translation energy modes are assumed in equilibrium at temperature T v . The species thermodynamic curve fits are 
still employed in this model by assuming that translation and rotation modes are fully engaged: 

_ ( m s + 2) R 
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B. Source Terms 
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where a sr and P sr are stoichiometric coefficients for species s as reactants and 


products, respectively, in reaction r . Reaction sets are taken from References 12 and 13. The reaction rate 
coefficients 12 ' 1 ' are functions of temperature: k f r = C r T 0 r exp(— E r /#T q ) and k b r = k f r (T) / K c r . The 

equilibrium constant, K c r , is derived from thermodynamic curve fit data for each species in the reaction. The 
effective temperature for a reaction, T q , is equal to T in the case of thermal equilibrium and exchange / shuffle 

reactions in thermal nonequilibrium. It is equal to T 0 ' 7 T v ° ' for dissociation reactions and equal to T v for electron 

impact ionization as an empirical approximation to the effect of thermal nonequilibrium in these processes. A 
minimum rate controlling temperature may be specified to preclude stiffness related convergence problems. These 
source terms and their Jacobians with respect to conserved variables are computed within the HEFSS kinetics 
module. The effect of turbulence models on kinetic rates 4 are not yet modeled. 
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Energy mode source terms: W = 


— n e S I S . The first term represents the relaxation 


between vibration-electronic modes and translation-rotation modes (combine terms 6 and 7 of Equation 16 in 

* 

Reference 11) and e v is the vibration-electronic energy for species s at thermal equilibrium conditions 
(T = T) . The second term represents the energy loss due to electron impact ionization. 11 The depletion of 
vibrational energy through dissociation uses the term e vs to model the (possibly) preferential vibration energy level 

of dissociating molecules. These source terms and their Jacobians with respect to conserved variables are computed 
within the HEFSS thermal relaxation module. 


C. Molecular Transport 

Transport properties for a perfect gas are calculated from Sutherland’s law and constant Prandtl number. Curve 
fits are used for equilibrium air 14 . Polynomial fits for collision cross sections of species pairs as a function of 
temperature 15 are used to compute transport properties for gas in chemical nonequilibrium. 11 The binary diffusion 
approximation is employed in Eq. 1 . The summation of the diffusion term in Eq. 1 over all species should equal 
zero. A mass weighted correction term as defined in Reference 16 enforces this zero sum condition. 
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D. Turbulence Models 

Linear eddy viscosity models assume that the shear stress angle and the mean flow angle are the same. Rumsey 
and Gatski 1 investigated the effect of K — Sand K — CO forms (both accommodated within HEFSS) of two- 
equation models on a multi-element airfoil. They demonstrated that models based on the K — 8 formulation were 
typically deficient in predicting adverse pressure gradient flows. Models based on the K — CO formulation 
provided a better prediction of adverse pressure gradient flows. Only the K — CO linear eddy viscosity model of 
Wilcox ls is tested here. 

The turbulent kinetic energy and specific dissipation rate are calculated by solving the following transport 
equations fully coupled to the species, momentum, and energy conservation: 
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where production term P and the coefficients a K , G c0 , y, 


P, fp. are derivable from those in Wilcox 18 . 


IV. Algorithm 

The reader should refer to FUN3D documentation 5 for a description of the baseline algorithm. The presentation 
that follows highlights algorithm modifications implemented for the hypersonic (high energy - FiEFSS) path 
through the software. 

Flux Reconstruction: Reconstruction refers to the algorithm used to derive the right and left states across a cell 
wall required in the flux definition in finite-volume, approximate Riemann solvers. Limiters refer to the algorithm 
used to constrain the reconstructed right and left states to prevent introduction of new extrema in the flowfield. 
Engagement of a limiter typically reduces the local order of accuracy of the reconstructed flux from 2 nd to 1 st order. 
Aftosmis et al. discuss the effects of various linear reconstruction techniques on unstructured meshes in Ref. 19. 
They note an increased degradation of solution quality associated with engagement of limiters on triangular (un- 
structured) meshes as compared to quadrilateral (structured) meshes. They also note that the use of directional 
limiting (local engagement of the limiter for each edge direction) as opposed to a global limiting (engagement of the 
limiter considering all directions emanating from a node) produced higher quality solutions. These observations are 
consistent with behavior observed in the FIEFSS benchmark simulations, to be discussed subsequently. 

The basic algorithm uses the divergence theorem (Green-Gauss) for viscous gradients and Roe’s method 20 for 
inviscid flux: 


f = 


1 
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[f L + f R 


# -1 |A|(dq - dq lim )] 


( 9 ) 


where 

dq = St( q L - q R ) (10) 

91 and 91 are the row and column eigenvector matrices of the flux Jacobian and A is the corresponding diagonal 
matrix of eigenvalues. The formulation of the Jacobian of f with respect to q is presented in Appendix A. Two 
options are offered here to reconstruct the left (L) and right (R) states. The first (Option 1) uses a least squares gra- 
dienL with a limiter function by Bartlr or Venkatakrishnan . The term dq lim in Eq. 9 is identically zero for 
Option 1. Option 1 is contained in the original FUN3D formulation. The second (Option 2) uses the nodal values for 
the left and right states. Option 2 achieves second-order spatial accuracy with the Symmetric Total Variation 
Diminishing 24 (STVD) scheme applied as a directional limiter (function of edge orientation) as follows: 
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and \j / is 0 if its argument is greater than 4 (fully first order), 1 if its argument is less than 2 (STVD limiting fully 

engaged), and linearly varying between these limits. The gradients Vq M at nodes (M = L or R) are computed using 
the identical, unlimited least squares formulation of primitive variables from Option 1 . The STVD limiting in Option 
2 and the Venkatakrishnan limiting in Option 1 are observed to diminish but not eliminate overshoots and 
undershoots in temperature and pressure in the vicinity of strong shocks. The pressure ratio weighting function, 4>, is 
used to smoothly switch to first-order in the vicinity of strong shocks and under-resolved expansions and to augment 
the limiting process in Option 2. Options 1 and 2 use identical formulation of the viscous terms. 

Eigenvalue Limiting: Elements, A |ct|n , of the diagonal eigenvalue matrix |A| in Eq. 9 are limited from 

below (entropy fix) 25 to enhance numerical dissipation across edges where a projected wave speed (Appendix A) 
approaches zero. 


^-ref - A, 0 ( a + U) 


Here, 0 < A, 0 < 1 for inviscid flows. Application of this limiter across the boundary layer adversely affects 

computed surface heating and skin friction. 26 It is easily scaled down in structured grid codes (LAURA) across the 
boundary layer. A requisite rescaling in the unstructured grid context is based on an evaluation of a local edge 
Reynolds number where the reference length is a function of radii of inscribed spheres contained in tets surrounding 
the edge. A loop over cells is executed to compute the inscribed radius for each cell. This value, 
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r ins node . The edge Reynolds number is defined Re edge = (par ins / |u) edge where edge values are averaged from 

the left and right nodal values. The value of X () is reduced by a factor min[l, (Re edge / ioo ) 7 8 J in viscous flows. 

The value of edge Reynolds number across a reasonably resolved boundary layer is expected to vary between 1 and 
100; consequently the eigenvalue limiting is effectively turned off in resolved boundary layers and shear layers. 

Grid Adaptation: A simple adaptation algorithm has been introduced to reduce cell size and accommodate high 
aspect ratio tets aligned with captured shocks. The intent is to reduce effects of “stair stepping” as non-aligned, 
isotropic grids produce “jagged” shocks, leading to erratic entropy fields following the streamlines that cross them. 
The algorithm is intended to serve as a precursor to adjoint-based adaptation that is available in the baseline 


FUN3D. A tensile force on each edge is defined F RL = 1 + K 


Feature based adaptation is 


(gRgL ) W2 

controlled by a user specified function g . The displacement at each node is defined by the sum of the tensile forces 
along each edge connected to the node. The undamped expression for the displacement of node L is defined 
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Displacements are damped by a factor 0.5 and limited by the value of 0.5r ins node before updating interior nodes. 

Updates at boundary surface nodes are constrained to move along the surface. (Grid points are not allowed to move 
on solid surfaces to prevent corruption of the original CAD surface which is not available as a reference in this 
simple approach.) The limitation by the inscribed radius is intended to prevent formation of negative volumes in 
high aspect ratio tets. The adaptation step is applied at a frequency of once per 20 to 50 global relaxation steps. 


V. Results and Discussion 


Inviscid Flow and Entropy Fix: Simulation of an 
erroneous carbuncle flow from the stagnation region of a 
blunt body in hypersonic flow with Roe schemes on 
structured grids can arise from a lack of natural 
dissipation in extensive regions where the eigenvalues 
are very small. The underlying scheme looks like a 
central difference formulation with associated odd-even 
decoupling. The problem is most evident when 
structured grids are closely aligned with the captured 
bow shock. The problem is overcome by limiting the 
minimum value an eigenvalue in Eq. 9 may assume, as 
specified in Eq. 12. Originally formulated to suppress 
rarefaction shocks, the eigenvalue limiting is often called 
an entropy fix. 

Initial tests of hypersonic flow over blunt bodies with 
unstructured grids did not show evidence of a carbuncle. 
It was surmised that natural dissipation in the scheme 
associated with edges that were not generally well 
aligned with the flow was sufficient to overcome this 
problem. However, close inspection of the stagnation 
region in an inviscid, hypersonic flow over an ellipsoid 
revealed irregularities still existed when Eq. 12 was not 
used to limit the eigenvalues. 

Fig. 1 shows surface pressure distributions in the 
stagnation region of an ellipsoid for a perfect gas, Mach 
6 flow at 30 degrees angle of attack. The simulations are 
identical except that the simulation displayed on the left 
employed almost no limiting while the simulation on the 
right used = 1 in Eq. 12. The limiting acts as a 

smoother in regions where it is engaged. It responds to 
the lack of numerical dissipation in the scheme in 
regions where eigenvalues are small compared to other 
relevant wave speeds (acoustic, convective). Clearly, 
some irregularities (non-elliptical shapes) persist in the 
contours even with the limiter. These are related to 
uneven shock capturing with a non-aligned grid 
(discussed in a later section). Though not shown here, 
solution convergence in leeside flows over bodies at 
high angle of attack is observed to benefit from the 
eigenvalue limiting process. In general, eigenvalue 
limiting is regarded as an essential aspect of the 
numerical formulation with Roe’s method for blunt 



Fig. 1: Effect of eigenvalue limiter on stagnation 
region pressure contours, over ellipsoid. 


HEFSS 


LAURA 



Fig. 2: Structured grid and equivalent uniformly 
biased unstructured grid in symmetry plane over 
cylinder for hypersonic test case. 
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body, hypersonic flow simulation on unstructured grids, 
even though the carbuncle phenomenon is not specifically 
observed. Restriction of the limiters in boundary layers 
and free shear layers to prevent corruption of computed 
heating levels is discussed below. 

Heating Benchmarks: Benchmark cases for 

hypersonic flow over a sphere and a cylinder at V M = 5000 
m/s, p x = 0.001 kg/m 3 , T x = 200 K (approximately Mach 
17) consistent with earlier LAURA benchmarks, have been 
intensively worked during the HEFSS code development 
process. The structured LAURA grid, adapted to align with 
the bow shock, is shown in the bottom part of Fig. 2. The 
unstructured grid for the HEFSS tests (top of Fig. 2) is 
formed directly from the structured grid. Both grids are 
three-dimensional (structured hexahedra (hexes), 
unstructured tetrahedra (tets)). An unstructured surface 
grid with ten spanwise cells is shown in Fig. 3. The grid is 
biased with diagonals running in the same direction. The 
biasing is intentionally retained to test simulation quality Fig. 3: Uniformly biased, unstructured surface 
under adverse conditions. Spanwise symmetry of the S r ' c * c * er ' vet * from structured grid with ten 
solution across the bias is easily monitored. Comparisons spanwise cells, 

to the structured grid LAURA results are used as 
benchmarks in the evaluation process. 

The single spanwise cell simulation results with Option 2 (STVD) are shown in Fig. 4. In this case there are no 
interior nodes in the spanwise direction. The HEFSS results for heating and pressure at surface nodes (indicated by 
symbol location) on the front and back plane are presented; the pressure symbols on the front and back plane exactly 
match, the heating levels exhibit a very small offset (less than half a symbol height). Good agreement with the 
LAURA reference is achieved and the near perfect over plotting of the front and back plane results indicates 
excellent spanwise constancy. The ten spanwise cell simulation results with Option 2 (STVD) are shown in Figs. 5 
and 6. In this case there are nine interior nodes on each spanwise line as shown in Fig. 3. The value at every surface 
node from the HEFSS simulation is plotted. The eleven spanwise nodes for pressure again show excellent overlap 
but the heating shows significant spanwise variation approaching the y=0 plane (Fig. 5) with the maximum spread in 
the stagnation region and solution quality improving away from the stagnation region. The average value of both the 
pressure and heating continue to be in good agreement with the LAURA results. The equivalent result using Option 
1 with the Venkatakrishnan limiter is shown in Fig. 7. In this case, pressure constancy is slightly degraded and more 





Fig. 4: Pressure and heating distribution for Fig. 5 . Heating contours in stagnation region of 

cylinder test case with 3D, uniformly biased, cylinder with 10 spanwise cells 

unstructured grid across one spanwise cell. 
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Fig. 6: Pressure and heating distribution for 
cylinder test case with 3D, uniformly biased, 
unstructured grid across ten spanwise cells with 
Ontion 2. 


Fig. 7: Pressure and heating distribution for 
cylinder test case with 3D, uniformly biased, 
unstructured grid across ten spanwise cells with 
Option 1. 


significant asymmetry in heat transfer is evident. Though not shown, the Option 1 results with the Barth limiter are 
of poorer quality and appear more like a first-order residt. In all cases, the symmetry of shock layer contour plots of 
pressure, temperature, and Mach number show consistent trends; the spanwise symmetry properties of the Option 2 
path are superior to all other tested options. 

These results demonstrate that solution quality for stagnation point heating behind strong shocks using tets all 
the way to the wall is poor using any of the reconstruction options discussed herein. It is unfair to lay all of the 
blame on the Green-Gauss formulation of viscous gradients across high aspect ratio cells when the results are so 
obviously sensitive to the higher order inviscid formulation. In any case, the freedom to utilize tets across the entire 
domain - even to the wall - provides the greatest flexibility for future grid adaptation. Note that the first-order, 
edge-based reconstruction in Eq. 9 only involves information at nodes R and L. A conjecture is offered that a truly 
multi-dimensional reconstruction of the flux at every 
face coidd overcome some of the simulation quality 
issues noted above. A multi-dimensional upwind (MDU) 
algorithm for high speed flows has recently been 
developed 27,28 ; it is not yet known how it would perform 
on the ten-spanwise-cell cylinder test problem for 
hypersonic, blunt-body heating. 

Reacting Gas Benchmark: The heating benchmark 

case described above has been repeated for a five species 
air model in chemical nonequilibrium and thermal 
equilibrium. Only the single spanwise cell grid is used. 

A fully catalytic wall boundary condition is 
approximated by specifying mole fractions at the wall 
equal to mole fractions in the cold freestream. The 
equivalent LAURA solution is again utilized as a 
benchmark. Excellent agreement between LAURA and 
HEFSS predictions for species mass fraction across the 
shock layer and boundary layer are observed in Fig. 8. 

Surface heating, including conduction and diffusion 
components, show excellent agreement (Fig. 9) outside 
of the stagnation region. The HEFSS prediction of 
stagnation point heating is slightly less than LAURA, 
consistent with results of Fig. 4. Shear stress is 
compared in Fig 10; it indicates an obvious sensitivity to 



Fig. 8: Species mass fraction distribution across the 
stagnation streamline. 
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Fig. 10: Shear stress distribution around cylinder 
at Mach 17. 



1.5 


0.5 


Fig. 9: Heating distribution around cylinder at 
Mach 17 for reacting air, fully catalytic wall. 

grid biasing in that asymmetry is obvious 
between the left and right side HEFSS results. 

The peak shear from HEFSS is greater than the 
LAURA result. Investigation of this issue is 
ongoing. 

Turbulent Flat Plate Benchmark: Turbulent 

flow of a perfect gas at Mach 6 and Reynolds 
number 26.38 million over a plat plate is 
simulated with HEFSS and VULCAN 
(benchmark). The wall to total temperature ratio 
is 0.6. Both codes used the same K-ro model 
without compressibility corrections. Transition 
is initiated at the leading edge of the plate. 

Comparisons of heating and skin friction 
coefficient are presented in Fig. 11. The skin 
friction coefficients are in excellent agreement 
except at the leading edge where boundary 
condition implementation differs between a 
node based (HEFSS) and cell center based 
(VULCAN) codes. The HEFSS heating is about 
10% lower than VULCAN over most of the 
plate. 

Grid Adaptation Test: The grid adaptation 
algorithm described in the previous section is 
tested on the inviscid ellipsoid (Fig. 1). The 
function g that controls tensile force F RL in Eq. 

(13) was keyed to pressure, velocity and temperature. An adaptation frequency of once per 20 relation steps was 
applied over approximately 5000 steps. A before and after view of the grid on the symmetry plane is shown in Fig. 
12. The grid clearly shows the extent of the captured bow shock. The clustering is accommodated with high aspect 
ratio cells. The corners of the domain are deformed because the algorithm to constrain movement on boundaries 
smoothes the surface normal at common boundaries. The cleaner shock capture results in significant improvement in 
the shape of pressure contour lines (Fig. 13) in the stagnation region. (Contrast Fig 13 with Fig. 1.) 


Fig. 11: Skin friction coefficient and heat transfer on flat 
plate at Mach 6 and Reynolds number 26.38 million. 
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Fig. 12: Original and adapted grids over ellipsoid in Fig. 13: Pressure contoui s over ellipsoid with 
symmetry plane. and without grid adaptation. 


Shuttle Orbiter: A simulation of hypersonic flow in thermochemical nonequilibrium (5-species air model) over 
the shuttle orbiter is presented to demonstrate present capabilities using only tets (1.1 million nodes) for 
discretization elements. The original grid was assembled very quickly (1 day) during the Columbia accident 
investigation using GridEx 29 software. The simulation was abandoned early on because it became clear that solution 
quality was not adequate with the existing algorithm. The simulation was revisited for this paper in order to 
document current capabilities on a reasonably complex problem, especially regarding need for grid adaptation. To 
repeat a theme developed earlier, we must demonstrate a capability to quickly generate an adequate grid and depend 
on the software to automatically adapt and refine it. In this case, the viscous near wall grid spacing was based on an 
existing structured grid solution from LAURA 30 at Mach 24. (STS-2, V,= 6920 m/s, p x = 5.75 I O ' kg/m 3 , T., = 
202 K, a* = 39.4 deg., and h = 72.4 km) The grid was marched out from the body in 40 layers. An essentially 
isotropic grid filled the volume between the viscous, near wall grid and the specified mesh spacing along the outer 
boundaries (farfield, exit, and symmetry plane). 




Fig. 14: Temperature contours in plane of 
symmetry over the nose in the original grid. 


Fig. 15: Temperature contours in plane of 
symmetry over the nose in the adapted grid. 
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Fig. 16: Temperature contours in plane of 
symmetry and exit plane behind the trailing edge 
in the original grid. 


Fig. 17: Temperature contours in plane of 
symmetry and exit plane behind the trailing edge 
in the adapted grid. 


The grid over the nose in Fig. 14 is clearly inadequate to cleanly capture the bow shock. However, sufficient 
resolution is provided to enable a sharper capture with high aspect ratio tets in Fig. 15. Further downstream, it is 
evident that the original grid (Fig. 16) has extremely poor resolution of the wing wake and leeside (no tail in this 
simulation) over the fuselage. The captured bow shock front is diffuse and wavy. The adapted grid (Fig. 17) again 
provides a cleaner shock capture and picks up an oblique shock over the deflected body flap (the back end of the 
shuttle is filled in over the body flap). However, the simple adaptation algorithm is unable to pull in more nodes to 
resolve vortex structure on the leeside or the wing wake. In fact, the leeside symmetry plane was plagued by an 
instability that is believed to be caused by very poor resolution of the separated vortex. A tendency for reverse flow 
at the exit plane between the OMS pods was also encountered; the solution would be reset at locations where this 
occurred as indicated by a local temperature falling outside a range of 0. IK to 30,000K. 

The ability to produce high aspect ratio cells aligned with the captured shock is a great advantage of this simple, 
spring-based adaptation strategy. However, the strategy as currently implemented cannot perform edge swapping 
nor can it provide local enrichment of nodes as obviously required on the wake and leeside. Isolated locations were 
observed where the grid appeared to be 
drawn to a point in the absence of any 
apparent driving gradients. This behavior is 
thought to be a remnant of the grid 
evolution process and restrictions to grid 
movement when the inscribed radius of any 
tet connected to a node gets small. (In the 
present case, no movement was allowed at 
any node where the inscribed radius was 
less than 0.001 inches.) Recent adaptation 
algorithms developed by Park 31 will be 
applied here in an effort to correct these 
deficiencies. 

Winds ide pressure contours from 
HEFSS on the adapted grid are compared 
to the LAURA solution in Fig. 18. 

Agreement is generally good. The 
unstructured grid solution included elevon 
gaps and body flap that were not included 
in the LAURA solution, thus accounting 
for differences at the trailing edge. The 


p, N/m 2 



Fig. 18: Windside pressure contours from unstructured 
HEFSS and structured LAURA. 
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Fig. 19: Leeside pressure contours from unstructured HEFSS 
and structured LAURA. 


nose region is in excellent agreement. The 
original, pre-adapted grid solution was very 
erratic looking - much like the behavior 
shown in Fig. 13 on the ellipsoid. 

The leeside pressures in Fig. 19 are in 
poor agreement, as woidd be expected given 
the poor leeside resolution discussed earlier. 

The compression over the canopy and vortex 
impingement on the OMS pod is greatly 
dissipated. 

Windside heating distributions are 
compared in Fig. 20. In this case, contour 
lines should not be expected to be mirror 
images because the surface boundary 
conditions are different. The HEFSS 
simulation used a constant wall temperature 
of 500 K and a fully catalytic wall boundary 
condition; it does not yet have the option for 
the finite catalytic, radiative equilibrium 
boundary condition used in LAURA. These 
differences do not have a strong influence on 
windside pressure in Fig. 18. The heating on 
the HEFSS simulation should be higher than 
LAURA as generally observed in Fig. 20 
with these boundary conditions. However, 
the HEFSS contour lines are more erratic in 
nature. Heating near the plane of symmetry 
shows a wake like behavior downstream of 
the stagnation region. The wake of the bow 
shock - wing shock interaction is more 
sharply defined in the HEFSS solution 
though contour lines remain somewhat 
jagged. The high heating around the elevon 
gaps is exaggerated because of surface grid is 
too coarse to resolve flow drawn into the gap. 

The results of Fig. 20 highlight the 
challenges remaining to bring HEFSS up to 
the standards required for vehicle design and 
analysis. Recall that we intentionally restrict 
all simulations to use only tets. This 
restriction arises from a belief that if semi- 
structured grid (hexes and prisms) are 
required in the near wall layers to achieve 
good / smooth heating then one is likely 
masking problems elsewhere in the flowfield 

- especially if separated flows or interactions off the surface are present. The assertion here is that if one uses semi- 
structured, high-aspect-ratio elements near the wall then one must also provide an algorithm to automatically 
generate and align equivalent elements elsewhere in the field as a function of flow physics. Such a grid adaptation 
algorithm may in fact be the best approach to this problem. The approach pursued here is to develop an algorithm 
that can provide a quality aerothermal simulation using only tets. The results in this paper highlight our best efforts 
to that end at this point. Several strategies adopted here (STVD limiting, entropy fix scaled to local metric of cell 
Reynolds number, high-aspect-ratio grid adaptation) have improved aerothermal simulation for hypersonic flows. 
The fact that the single spanwise cell cylinder results yield good aerothermal simulations lead us to believe that the 
problem is not so much that the tets cannot support a good solution but that current algorithms admit discretization 
errors that manifest themselves more severely in the randomly oriented unstructured grid context than in the 



Fig. 20: Windside heating contours from unstructured 
HEFSS and structured LAURA. 
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structured grid context. We speculate that the quasi-one-dimensional, edge oriented, reconstruction used here is the 
predominant cause of the problem and future work will be directed at this concern. 

VI. Concluding Remarks 

The synthesis of physical models for gas chemistry and turbulence from the structured grid codes LAURA and 
VULCAN into the unstructured grid code FUN3D is described. This suite of modules for hypersonic application is 
called the High Energy Flow Solver Synthesis (HEFSS). A Symmetric, Total Variation Diminishing (STVD) 
algorithm and entropy fix (eigenvalue limiter) keyed to local cell Reynolds number are introduced to improve 
solution quality for hypersonic applications focused on aeroheating. A simple grid-adaptation procedure is 
incorporated within the flow solver to reduce effects of jagged captured shocks on shock layer profiles. Simulations 
of flow over an ellipsoid (perfect gas, inviscid), Shuttle Orbiter (viscous, chemical nonequilibrium) and comparisons 
to the structured grid solvers LAURA (cylinder, Shuttle Orbiter) and VULCAN (turbulent flow on a flat plate) are 
presented to show current capabilities. 

Early in the testing phase of HEFSS it became apparent that stagnation region heating was poorly predicted. The 
poor predictions are much more a result of inviscid algorithm deficiencies than shortcomings associated with 
exclusive use of tetrahedral elements across the shock and boundary layer. Two sets of numerical results support this 
assertion. First, single-span wise-cell cylinder tests on a strongly biased, unstructured grid for both single-species and 
five-species flow in thermochemical nonequilibrium compared well with heating benchmarks from a structured grid 
code. In contrast, the ten-spanwise-cell cylinder tests (same grid as the single-spanwise-cell case repeated across the 
span) showed a marked degradation in heating contours in the stagnation region, apparently arising from spanwise 
flow admitted with the additional degrees of freedom. If grid were the primary culprit, one would not expect a good 
solution in the single-spanwise-cell tests. Second, numerical tests on the inviscid flux formulation (eigenvalue 
limiting, reconstruction options) had significant effects on the stagnation region heating with little effect outside the 
stagnation region. These results indicate that a primary cause of heating degradation arises from how the inviscid 
flow gets processed crossing the shock and the shock layer to the boundary-layer edge; and not from poor 
formulation of viscous terms across the boundary layer itself. Indeed, inviscid flow tests on an ellipsoid showed 
pressure contour irregularities (occurring only in the stagnation region) were fixed by application of an entropy fix 
and better grid alignment. If pressure (and entropy) at the boundary layer edge is not correct, there is no chance to 
get good prediction of heating. 

A Symmetric Total Variation Diminishing (STVD) algorithm with cell Reynolds number scaled eigenvalue 
limiting produced the best quality heating in the stagnation region of several approaches tested. However, all 
inviscid formulations tested utilized a quasi-one-dimensional reconstruction and none of them produced 
exceptionally good predictions of heating in the stagnation region. The STVD approach as applied herein will admit 
new maxima and minima across a captured shock and requires application of an auxiliary pressure switch for 
stability. Research continues on defining an inviscid reconstruction algorithm that will yield consistently good 
hypersonic stagnation heating using only tetrahedral elements across the entire shock layer. 
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Appendix A: Characteristic Variables and Jacobians 

Formation of characteristic variables and inviscid flux Jacobians closely follows that in LAURA 11 . However, 
the reconstruction variable set uses primitive variables rather than conservative variables. The Jacobian of the 
inviscid flux vector with respect to conserved variables is expressed below in Eq. A1 as a function of its 
eigenvalues. Note that subscripts s and r refer to species equations, subscripts i and j refer to momentum equations, 
and subscripts v and v refer to supplemental energy mode equations. A single, total energy equation is also included. 
The two equation turbulence models are not explicitly included here. Their form looks very much like the species or 
energy mode rows. Additional terms arise in the momentum and energy equations if turbulent pressure in the 
Reynolds stress tensor is treated as part of the inviscid flux contribution to pressure (rather than a source term which 
is not involved in the upwind discretization). 
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The eigenvalues of the Jacobian A are: A,j = U = U^ and X 23 = U + a . The functions of eigenvalues for 
the Jacobians A and |A| are now defined: 




jAi 


A 




A-i = 

IN 

for 

|A| 


(A2) 



| 0 


A 




A = 

{(2 Aj| - |A, 2 - |A, 3 ) /(2a 2 ) 

for 

|A| 


(A3) 



f 1 


A 




A = 

|(|A 2 | - A 3 ) /(2a) 

for 

|A| 


(A4) 

The elements of dq in Eq. 10 are defined: 





d q s = 
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where we 

have used dp = P(dpE 

: - UjdpUj) + cp u dpe u + y r dp r 

(again neglecting optional inclusion of 


turbulent pressure term P t pKas part of the upwind discretization). 
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